function ctrl = GuidanceStage2(statec, ctrlp)
global R0 m alphaMax alphaMin alphaLD S kh1 d2r
rad  = statec(1);
V = statec(4);
fpa = statec(5);
alphap = ctrlp(1);
[~,sonic,~,rho] = atmoscoesa(rad - R0);
Ma = V/sonic;
[CL, ~] = getCLCD(alphaLD, V/sonic);
[~, CD] = getCLCD(alphap, V/sonic);
q = 0.5*rho*V^2;
L = q*S*CL;
D = q*S*CD;
dhc = - (D/m) * 2 * 7110 / V;
% dhc = -2*g*Hs/LtoD/cos(bank)/V;
dh = V * sin(fpa);
Ld = L/m + kh1 * (dhc - dh);
fun = @(alpha)abs(Ld - ...
      q*S/m*(-0.21 + 0.005059*Ma + 0.0497*alpha/d2r + 0.0004244*(alpha/d2r).^2));
alphad = fminbnd(fun, alphaMin, alphaMax);
ctrl(1) = alphad;
ctrl(2) = 0;
ctrl(3) = 1;
end